%%  SelectSolver(cfg,exvars,depvars,indepvars,r,t,calibrationnames,startingvalues,varargin)
% Used in model calibration, simplifies the calling of solvers given their
% inconsistent notations.
% Calibrationnames: the variables to be optimized over
% Startingvalues:   their starting values
% varargin{1}: the Solver to be used (optional, default is gui setting)
% varargin{2}: the boundary conditions (B x <= A), with a the lowerboundaries
% (optional)

function [cfg,depvars,indepvars] = SelectSolver(cfg,exvars,depvars,indepvars,r,t,calibrationnames,startingvalues,varargin)
global globals
%% Determine the requested varargin options
for n = 1:length(calibrationnames)
    eval(['CalibrationLogical{n} = ones([1 size(' calibrationnames{n} ',3)]);' ]);
end

if ~exist('varargin','var') && isempty(varargin)
    Solver=cfg.Solver;
    lowerboundaries = zeros(length(startingvalues),1);
    upperboundaries = ones(length(startingvalues),1)*inf;
elseif length(varargin)==1
    Solver=varargin{1};
    lowerboundaries = zeros(length(startingvalues),1);
    upperboundaries = ones(length(startingvalues),1)*inf;

elseif length(varargin) ==2
    Solver=varargin{1};
    lowerboundaries = permute(varargin{2}, [3 1 2]);
    upperboundaries = ones(length(startingvalues),1)*inf;
elseif length(varargin) == 3
    Solver=varargin{1};
    lowerboundaries = permute(varargin{2}, [3 1 2]);
    upperboundaries = permute(varargin{3}, [3 1 2]);
elseif length(varargin) == 4
    Solver=varargin{1};
    lowerboundaries = permute(varargin{2}, [3 1 2]);
    upperboundaries = permute(varargin{3}, [3 1 2]);
    CalibrationLogical = varargin{4};
end
tic;

cfg.CalibrationArgs{r,1}= calibrationnames;
calibrationvector       = createvector(depvars,indepvars,exvars,calibrationnames,r,CalibrationLogical); % creates a vector with variable x and calibrationvariable

%% Save the variables to global for inspection of the process later
globals.Error             = [];  
for n=1:length(cfg.CalibrationArgs{r,1}); eval(['globals.CalibrationFactors.' (cfg.CalibrationArgs{r,1}{n}) '= [];']); end
if strcmp(Solver,'None')==1
    disp('Skipping optimization');
else
switch Solver
     case 'fmincon'
        %% Fmincon
         eval(['[x,depvars.ObjectiveFunction(r)] '...
         '= fmincon(@(x) ' cfg.Model '(cfg,exvars,depvars,indepvars,r,t,' calibrationvector ')'...
         ',startingvalues, -eye(length(startingvalues),length(startingvalues))'... 
         ',-lowerboundaries,[],[],[],[],[],cfg.MinOptions);']);

    case 'fminsearch'
        %% Fminsearch
        eval(['[x,depvars.ObjectiveFunction(r)] '...
        '= fminsearch(@(x) ' cfg.Model '(cfg,exvars,depvars,indepvars,r,t,' calibrationvector ')'...
        ',startingvalues,cfg.MinOptions);']);
    
    case 'fminbnd'
        %% Fminbnd
        eval(['[x,depvars.ObjectiveFunction(r)] '...
        '= fminsearchbnc(@(x) ' cfg.Model '(cfg,exvars,depvars,indepvars,r,t,' calibrationvector ')'...
        ',lowerboundaries,[],cfg.MinOptions);']);
    
    case 'fminsearchbnd'
        %% Fminsearchbnd
        eval(['[x,depvars.ObjectiveFunction(r)] '...
        '= fminsearchbnd(@(x) ' cfg.Model '(cfg,exvars,depvars,indepvars,r,t,' calibrationvector ')'...
        ',startingvalues,lowerboundaries,upperboundaries,cfg.MinOptions);']);
        
    case 'lsqnonlin'
        %% Lsqnonlin
        eval(['[x,bs1,depvars.ObjectiveFunction(r)]  = lsqnonlin(@(x) ' cfg.Model '(cfg,exvars,depvars,indepvars,r,t,' calibrationvector ')'...
            ',startingvalues);']);

    case 'ga'
        %% Genetic Algorithm
        eval(['[x,depvars.ObjectiveFunction(r)] = ga(@(x) ' cfg.Model '(cfg,exvars,depvars,indepvars,r,t,' calibrationvector ')'...
            ',length(startingvalues),-eye(length(startingvalues),length(startingvalues))'...
            ',-lowerboundaries,[],[],lowerboundaries,upperboundaries);']);
        
    case 'MultiObj Genetic Algorithm'
        %%
        X = gamultiobj(FITNESSFCN,NVARS);

    otherwise
    disp('Invalid Solver');

end

%% Convert the output into meaningful variables
disp(['Optimization completed in ' num2str(toc) 's, with ' Solver ', Objective Function = ' num2str(depvars.ObjectiveFunction(r))]);
temp.decomposedvarargin = textscan( calibrationvector,'%s','delimiter',',');
for n = 2:2:size(temp.decomposedvarargin{1},1)
    eval( [eval(temp.decomposedvarargin{1}{n-1}) '(r,1,:)=' temp.decomposedvarargin{1}{n} ';']);   
end
end

end


function vectorout = createvector(depvars,indepvars,exvars,calibrationnames,r,CalibrationLogical)
vectorout = [''];
m         = 1;
for n=1:length(calibrationnames)        
        vectorout(end+1:end+length(calibrationnames{n})+2) = ['''' calibrationnames{n} ''''];
        if eval(['size(' calibrationnames{n} ',3)']) == 1
            vectorout(end+1:end+5+ceil(log10(m+0.001))) = [',x(' num2str(m) '),']; m=m+1;
        elseif eval(['size(' calibrationnames{n} ',3)']) > 1
            vectorout(end+1:end+2) = ',[';
            for i = 1:eval(['size(' calibrationnames{n} ',3)'])
                if indepvars.UseModes(r,1,i) && CalibrationLogical{n}(i)==1
                    vectorout(end+1:end+4+ceil(log10(m+0.001))) = ['x(' num2str(m) ') ']; m=m+1;
                else
                    vectorout(end+1:end+2) = '0 ';
                end
            end
            vectorout(end:end+1)='],';
        else
            error('Invalid third dimension of matrix');
        end
end
vectorout = vectorout(1:end-1);
end